package odesimulation;

import java.awt.GridLayout;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.HashSet;
import java.util.Iterator;

import javax.swing.JFrame;
import javax.swing.JPanel;

import odesimulation.crnParser.ReactionsetContext;

import org.antlr.v4.runtime.ANTLRInputStream;
import org.antlr.v4.runtime.CommonTokenStream;
import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.nonstiff.ClassicalRungeKuttaIntegrator;
import org.apache.commons.math3.ode.nonstiff.GraggBulirschStoerIntegrator;

//import picoevo.app.chemicalreactions.Global_CR;
//import picoevo.app.chemicalreactions.Individual_CR;

import fitness.Constants;
import fitness.ODESineFitness;

import utils.MyStepHandler;
import utils.PlotExpData;
import utils.TargetBifurcationDiagram;

public class OdeImplicitSolver_Nicolas implements FirstOrderDifferentialEquations{

	public static double numberOfSteps = 1000;
	
	public static boolean debug = true;
	/**
	 * Like OdeSolver, but instead of using an explicit ode system, we only use
	 * an ArrayList of symbols and an ArrayList of Reactions on those.
	 */
	
	ArrayList<String> symbols;
	HashMap<String,Enzyme> enzymes; //map an enzyme to the list of stuff saturating it
	ArrayList<Reaction> reactions;
	
	static TargetBifurcationDiagram target = null;
	
	public OdeImplicitSolver_Nicolas(ArrayList<String> symbols, HashMap<String,Enzyme> enzymes, ArrayList<Reaction> reactions){
		this.symbols = symbols;
		this.enzymes = enzymes;
		this.reactions = reactions;
	}
	
	public static void setTarget(TargetBifurcationDiagram newTarget){
		target = newTarget;
	}
	
	public static void setDefaultTarget(){
		
		//target = TargetBifurcationDiagram.japaneseFlag();
		target = TargetBifurcationDiagram.twoBandsLR();
	}
	
	
	
	@Override
	public void computeDerivatives(double t, double[] y, double[] ydot) {
		
		Arrays.fill(ydot, 0); //just to be safe, put everything to 0
		HashMap<String, Double> conc = new HashMap<String, Double>();
		for(int i=0; i<y.length; i++){
			conc.put(symbols.get(i), y[i]);
		}
		//Update enzymes
		HashMap<String, HashMap<String, Double>> enz = new HashMap<String, HashMap<String, Double>>();
		for(String e :this.enzymes.keySet()){
			double partialSaturationValue = 1.0;//Still has to be devided in the final value by the relevant Km
			Enzyme actualEnzyme = this.enzymes.get(e);
			HashMap<String, Double> actualKis = new HashMap<String, Double>();
			for(String sat : actualEnzyme.Kms.keySet()){
				partialSaturationValue += conc.get(sat)/actualEnzyme.Kms.get(sat);
			}
			for(String sat : actualEnzyme.Kms.keySet()){
				actualKis.put(sat, actualEnzyme.kis.get(sat)/partialSaturationValue);
			}
			enz.put(e, actualKis);
		}
		
		for(Reaction r : reactions){
			HashSet<String> relevantSpecies = new HashSet(r.products); //so we remove double
			relevantSpecies.addAll(r.reactants);
			Iterator<String> it = relevantSpecies.iterator();
			while(it.hasNext()){
				String s = it.next();
				if(r.getClass() != EnzymaticReaction.class){
				ydot[symbols.indexOf(s)] += r.computeFlow(s, conc);
				} else {
					String enzyme = ((EnzymaticReaction) r).enzyme;
					Enzyme actualEnzyme = this.enzymes.get(enzyme);
					String origin = r.reactants.get(0);
					double RelevantKm = actualEnzyme.Kms.get(origin);
					double RelevantKi = enz.get(enzyme).get(origin);
					
					ydot[symbols.indexOf(s)] += ((EnzymaticReaction) r).computeFlow(s, RelevantKi/RelevantKm, conc);
				}
			}
		}
	}

	@Override
	public int getDimension() {
		return symbols.size();
	}

	public static ArrayList<Reaction> generateReactions(ReactionsetContext context){
		ArrayList<Reaction> reaclist = new ArrayList<Reaction>();
		
		return reaclist;
	}
	
	/**
	 * Solve a specific set of differential equations
	 * @param initialConditions Initial concentrations for all substrats
	 * @param symbols Names of substarts
	 * @param enzymes Enzymes present in the system. Link between their name and their physical representation
	 * @param reactions All the reactions present in the system.
	 * @return
	 */
	public static double[][] solveEqus(double[] initialConditions, ArrayList<String> symbols,  HashMap<String,Enzyme> enzymes, ArrayList<Reaction> reactions){
		
		OdeImplicitSolver_Nicolas odes = new OdeImplicitSolver_Nicolas(symbols, enzymes, reactions);
		GraggBulirschStoerIntegrator myIntegrator = new GraggBulirschStoerIntegrator(
				1e-13, 1, 1e-6, 1e-6);
		MyStepHandler handler = new MyStepHandler(odes.getDimension());
		myIntegrator.addStepHandler(handler);
		myIntegrator.addEventHandler(new utils.CalcErrEventHandler(), 2000, 0.000000001, 100000);
		myIntegrator.addEventHandler(new utils.MyODEEventHandler(), 10, 1, 100);
		myIntegrator.integrate(odes, 0, initialConditions, OdeImplicitSolver_Nicolas.numberOfSteps ,
				initialConditions);
		return handler.getTimeSerie();
	}
	
// =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

	static public void testDefault ( )	 // Example from Nath - 2014-01-09 
	{
	
		//Lexing the equations
		double k1 = 0.1;
		double maxk1 = 4.0;
		int sizex = 5;
		int sizey = 5;
		double mink4 = 0.0;
		double maxk4 = 8.0;
		double stepx = (maxk4-mink4)/((double)sizex);
		double stepy = (maxk1-k1)/((double)sizey);
	
		String exp = "";
		JFrame frame = new JFrame();
		JPanel content = new JPanel(new GridLayout(0,sizex));
		double time = System.currentTimeMillis();
		while(k1 < maxk1){
			double k4 = mink4;
			while(k4 < maxk4){
	
				crnLexer myLex = new crnLexer(new ANTLRInputStream(
					"# individual \"99\"\n"+
					"enzymes: E1 E2 E3"+
					"products: X Y Z"+
					"K_m[3][3] : { { 35.35651646265842 1.0 18.074538782097136 }{ "+(1.0/(k1*8.7e-2))+" 0.19639942990693068 1118.6216070905916 }{ 49.06140644009669 7.5529597327930205 4.656470904572082 } }"+
					"k_i[3][3] : { { 5392.936282100518 0.39 3.3851197232870134 }{ "+(1.0/8.7e-2)+" 6.127154985599243 7.02027575785981 }{ 275.49850207203826 16.44430374780971 2.3695397213682563 } }"+
						
					/**/
					// Nath default
					"X + E2 ==> E2 + X+X "+
					"X + Y <==> Y+Y [ k+ = 1.0 ; k- = 0.0 ]"+
					"X <==>   [ k+ = "+ 0.39*k4+" ; k- = 0.000 ]"+      // [n] Why modify K+? can we assume k+=1.0 ?  
					"Y + E1 ==>  E1"+
					" <==> X [ k+ = 0.0000001 ; k- = 0.0 ]"
					/**/
	
					/*
					// [Lotka1920] (Reference Oscillator )
					"X <==> X + X [ k+ = 1.0 ; k- = 0.0 ]"+
					"Y + X <==> Y + Y [ k+ = 1.0 ; k- = 0.0 ]"+
					"Y <==> [ k+ = 1.0 ; k- = 0.0 ]"
					/**/
	
					/*
					// [Lotka1920] + small value of k- = from k-=0.001 to k-=0.0001: less and less damped (but still damped)  ; k-=0.00001: not damped anymore
					// Changing the k+ value from 1.0 to 10.0 does not change anything. 
					"X <==> X + X [ k+ = 1.0 ; k- = 0.001 ]"+
					"Y + X <==> Y + Y [ k+ = 1.0 ; k- = 0.001 ]"+
					"Y <==> [ k+ = 1.0 ; k- = 0.001 ]"
					/**/
	
					/*
					// [Lotka1920] + minor change of k+ (1.0, 5.0, 10.0) value for first reaction = all oscillate, though with different shapes
					"X <==> X + X [ k+ = 10.0 ; k- = 0.0 ]"+
					"Y + X <==> Y + Y [ k+ = 1.0 ; k- = 0.0 ]"+
					"Y <==> [ k+ = 1.0 ; k- = 0. ]"
					/**/
	
					/*
					// [Lotka1920] + minor change of k- value for last reaction = damped oscillation 
					"X <==> X + X [ k+ = 1.0 ; k- = 0.0 ]"+
					"Y + X <==> Y + Y [ k+ = 1.0 ; k- = 0.0 ]"+
					"Y <==> [ k+ = 1.0 ; k- = 0.1 ]"
					/**/
	
					/*
					// [Lotka1920] + minor change of k+ (1.0, 5.0, 10.0) value for first reaction = all oscillate, though with different shapes
					"X <==> X + X [ k+ = 10.0 ; k- = 0.0 ]"+
					"Y + X <==> Y + Y [ k+ = 1.0 ; k- = 0.0 ]"+
					"Y <==> [ k+ = 1.0 ; k- = 0. ]"
					/**/
	
					/*
					// [Lotka1920] + adding E1 = minor change of the oscillation behavior.
					"X <==> X + X [ k+ = 1.0 ; k- = 0.0 ]"+
					"Y + X <==> Y + Y [ k+ = 1.0 ; k- = 0.0 ]"+
					"E1 + Y ==> E1"
					/**/
	
					/*
					// [Lotka1920] + adding E1 + adding Z = X,Y,Z (micro-)oscillate, Z goes up
					"X <==> X + X [ k+ = 1.0 ; k- = 0.0 ]"+
					"Y + X <==> Y + Y [ k+ = 1.0 ; k- = 0.0 ]"+
					"E1 + Y ==> E1 + Z"
					/**/
						));
	
	
	
				//Parsing the equations
				crnParser myPars = new crnParser(new CommonTokenStream(myLex));
				//Parsing result
				//Contains "symbols": X,Y,Z... list of all substrats.
				//		   "enzymes": E1... list of all enzymes.
				//         "reactions": list of all reactions.
				ReactionsetContext context = myPars.reactionset();
	
				//initial concentrations for substrats: 10 nM all around!
				double[] init = new double[context.symbols.size()];
				Arrays.fill(init, 5);
	
				double[][] result = solveEqus(init,context.symbols,context.enzymes,context.reactions);
				//System.out.println(k1+" "+k4);
				fitness.OscillatorBasedResult fit = ODESineFitness.evaluateFitness(result, context.symbols);
				exp += (fit.finalResult()>1?" "+fit.average+"\t":" -\t");
				String[] seqNames = new String[context.symbols.size()];
				seqNames =context.symbols.toArray(seqNames);
	
				PlotExpData demo = new PlotExpData("",
						result, seqNames, "", "");
	
				content.add(demo);
	
				k4 += stepx;
			}
			//System.out.println(exp);
			exp += "\n";
			k1+= stepy;
		}
		double time2 = System.currentTimeMillis();
		System.out.println(exp);
		System.out.println("\n"+(time2-time));
		frame.add(content);
		frame.pack();
		frame.setVisible(true);
		
	}
	
	
// =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
/*
	//Fitness to match flat individuals
	static public double getFitnessValue_1 ( Individual_CR individual ) // HACK: uses Individual_CR object from PicoEvoLite. Using a string representation would be better. (but must be parsed to change specific values)
	{

		Individual_CR indClone = (Individual_CR) individual.clone(); // build a clone that we can toy with.

		
		if ( debug )
			System.out.println(individual.getDescription());
		
		// [n]: IMPORTANT NOTE: I (nicolas) assume that only template initial concentration, K_m (and possibly k_i) can be modified. k+ and k- are *not* modified.
		
		//Lexing the equations
		double k1 = 0.1;
		double maxk1 = 4.0;
		int sizex = 5;
		int sizey = 5;
		double mink4 = 0.0;
		double maxk4 = 8.0;
		double stepx = (maxk4-mink4)/((double)sizex);
		double stepy = (maxk1-k1)/((double)sizey);
	
		double[][] matrixKm = indClone.getCopyOfMatrixKm(); // get Km values 
		double[][] matrixKi = indClone.getCopyOfMatrixKi(); // get Ki values
		
		// compute phase diagram
		
		double[][] bifurcationDiagram = new double[ (int) (( maxk4 - mink4 ) / stepx) ][ (int)(( maxk1 - k1 ) / stepy) ];
		int i = 0;
		int j = 0;
			
		while(k1 < maxk1){
			
			double k4 = mink4;
			i=0;			
			
			while(k4 < maxk4){
				
				matrixKm[1][0] = 1.0/(k1*8.7e-2); // overwrite
				matrixKi[1][0] = 1.0/8.7e-2; // overwrite
				
				indClone.setMatrixKm(matrixKm);
				indClone.setMatrixKi(matrixKi);
				
				crnLexer myLex = new crnLexer(new ANTLRInputStream( indClone.getDescription() ) );
	
				//Parsing the equations
				crnParser myPars = new crnParser(new CommonTokenStream(myLex));
				//Parsing result
				//Contains "symbols": X,Y,Z... list of all substrats.
				//		   "enzymes": E1... list of all enzymes.
				//         "reactions": list of all reactions.
				ReactionsetContext context = myPars.reactionset();
	
				//initial concentrations for substrats: 10 nM all around!
				double[] init = new double[context.symbols.size()];
				Arrays.fill(init, 5);
	
				double[][] result = solveEqus(init,context.symbols,context.enzymes,context.reactions);
				
				fitness.OscillatorBasedResult fit = ODESineFitness.evaluateFitness(result, context.symbols);
				
				bifurcationDiagram[i][j] = fit.finalResult();//>1?fit.average:0;
				
				k4 += stepx;
				i++;
			}
			
			k1 += stepy;
			j++;
		}

		// compute fitness (compare phase diagram with target phase diagram

		// TODO: toy fitness, target is flat everywhere
		double fitnessValue = 0; //Double.MAX_VALUE;
		for ( j = 0 ; j != bifurcationDiagram[0].length ; j++ )
		{
			for ( i = 0 ; i != bifurcationDiagram.length ; i++ )
				fitnessValue += bifurcationDiagram[i][j];
		}
		
		// debug
		if ( debug )
		{
			System.out.println("Bifurcation diagram: ");
			for ( j = 0 ; j != bifurcationDiagram[0].length ; j++ )
			{
				System.out.print(" ");
				for ( i = 0 ; i != bifurcationDiagram.length ; i++ )
					System.out.print(bifurcationDiagram[i][j]+"\t");
				System.out.println();
			}
			System.out.println("Fitness value: " + fitnessValue);
		}
		
		return fitnessValue;
	}
	*/
	
	//Fitness for normal evolution
	static public double getFitnessValue_2 ( String s) 
	{
		
		double[][] bifurcationDiagram = getDiagram(s, 5, 5);

		// compute fitness (compare phase diagram with target phase diagram

		// TODO: toy fitness, target is flat everywhere
		double fitnessValue = 0; //Double.MAX_VALUE;
		for ( int j = 0 ; j != bifurcationDiagram[0].length ; j++ )
		{
			for ( int i = 0 ; i != bifurcationDiagram.length ; i++ ){
				fitnessValue+= bifurcationDiagram[i][j];
			}
		}
		
		// debug
		if ( debug )
		{
			System.out.println("Bifurcation diagram: ");
			for (int j = 0 ; j != bifurcationDiagram[0].length ; j++ )
			{
				System.out.print(" ");
				for ( int i = 0 ; i != bifurcationDiagram.length ; i++ )
					System.out.print(bifurcationDiagram[i][j]+"\t");
				System.out.println();
			}
			System.out.println("Fitness value: " + fitnessValue +"\n");
		}
		
		return fitnessValue;
	}
	
	//Fitness for normal evolution
	static public double getFitnessValue_3 (String s, int res) 
	{
		
		int sizex = (int) Math.pow(2, TargetBifurcationDiagram.maxRes - res);
		int sizey = (int) Math.pow(2, TargetBifurcationDiagram.maxRes - res);
		
		if(target == null){
			setDefaultTarget();
		}
		
		double[][] bifurcationDiagram = getDiagram(s, sizex, sizey);
		double[][] resTarget = target.getDiagram(res);
		
		// TODO: toy fitness, target is flat everywhere
		double fitnessValue = 1; //Double.MAX_VALUE;
		for ( int j = 0 ; j != bifurcationDiagram[0].length ; j++ )
		{
			for ( int i = 0 ; i != bifurcationDiagram.length ; i++ ){
				double diff = resTarget[i][j] - bifurcationDiagram[i][j];
				fitnessValue+= diff*diff;
			}
		}
		
		fitnessValue = 1.0/fitnessValue; //We want to maximize
		
		// debug
		if ( debug )
		{
			System.out.println("Bifurcation diagram: ");
			for ( int j = 0 ; j != bifurcationDiagram[0].length ; j++ )
			{
				System.out.print(" ");
				for ( int i = 0 ; i != bifurcationDiagram.length ; i++ )
					System.out.print(bifurcationDiagram[i][j]+"\t");
				System.out.println();
			}
			System.out.println("Fitness value: " + fitnessValue +"\n");
		}
		
		return fitnessValue;
	}
	
// =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

	private static double[][] getDiagram(String s, int sizeX, int sizeY){
		
		int sizex = sizeX;
		int sizey = sizeY;
		
		//Lexing the equations
		crnLexer myLex = new crnLexer(new ANTLRInputStream(s));
		//Parsing the equations
		crnParser myPars = new crnParser(new CommonTokenStream(myLex));
		//Parsing result
		//Contains "symbols": X,Y,Z... list of all substrats.
		//		   "enzymes": E1... list of all enzymes.
		//         "reactions": list of all reactions.
		ReactionsetContext context = myPars.reactionset();
		
		//initial concentrations for substrats: 10 nM all around!
		// 10 nmol/L ?
		double[] init = new double[context.symbols.size()];
		
		
		if ( debug )
			System.out.println(s);
		
		// [n]: IMPORTANT NOTE: I (nicolas) assume that only template initial concentration, K_m (and possibly k_i) can be modified. k+ and k- are *not* modified.
		
		
		double k1 = Constants.range1min;
		double k2 = Constants.range2min;
		double stepx = (Constants.range2max-Constants.range2min)/((double) sizex);
		double stepy = (Constants.range1max-Constants.range1min)/((double) sizey);
		
		// compute phase diagram
		
		double[][] bifurcationDiagram = new double[ (int) (( Constants.range2max - Constants.range2min ) / stepx) ][ (int)(( Constants.range1max - Constants.range1min ) / stepy) ];
		int i = 0;
		int j = 0;
			
		while(k1 < Constants.range1max){
			
			k2 = Constants.range2min;
			i=0;			
			
			while(k2 < Constants.range2max){
				
				//overwrite params
				Arrays.fill(init, 5);
				setParams(context, init, k1, k2);
	
				double[][] result = solveEqus(init,context.symbols,context.enzymes,context.reactions);
				
				fitness.OscillatorBasedResult fit = ODESineFitness.evaluateFitness(result, context.symbols);
				
				bifurcationDiagram[i][j] = fit.finalResult();//>1?fit.average:0;  /Ftagn
				
				k2 += stepx;
				i++;
			}
			
			k1 += stepy;
			j++;
		}
		
		return bifurcationDiagram;
	}
	
	private static void setParams(ReactionsetContext context, double[] init, double k1,
			double k2) {
		String[] params = Constants.currentTarget1.split(" ");
		String enz, species="";
		Enzyme trueEnz;
		
		boolean debug = false;
		
		switch(Constants.typeCurrentTarget1){
		case KM:
			enz = params[0];
			if(params.length >= 2){
			species = params[1];
			}
			trueEnz = context.enzymes.get(enz);
			if(debug)
				System.out.println("Enz: "+enz+" species: "+species+" trueEnz: "+trueEnz);
			if(trueEnz != null && (species.equals("") || trueEnz.Kms.get(species)!=null)){
				if(debug)
					System.out.println("Done");
				if(params.length>=2){
					trueEnz.Kms.put(species, k1);
				} else { //we actually update everyone, with a multiplying factor
					for(String s : trueEnz.Kms.keySet()){
						//find the original value we stored
						double defaultValue = context.Kmvalues.get(context.enzymenames.indexOf(enz)).get(context.symbolList.indexOf(s));
						//multiply by the rangeValue
						trueEnz.Kms.put(s, k1*defaultValue);
					}
				}
			}
			break;
		case KI:
			enz = params[0];
			if(params.length >= 2){
			species = params[1];
			}
			trueEnz = context.enzymes.get(enz);
			if(debug)
				System.out.println("Enz: "+enz+" species: "+species+" trueEnz: "+trueEnz);
			if(trueEnz != null && (species.equals("") || trueEnz.kis.get(species)!=null)){
				if(debug)
					System.out.println("Done");
				if(params.length>=2){
					trueEnz.kis.put(species, k1);
				} else { //we actually update everyone, with a multiplying factor
					for(String s : trueEnz.Kms.keySet()){
						//find the original value we stored
						double defaultValue = context.Kmvalues.get(context.enzymenames.indexOf(enz)).get(context.symbolList.indexOf(s));
						//multiply by the rangeValue
						trueEnz.kis.put(s, k1*defaultValue);
					}
				}
			}
			break;
		case CONC:
			species = params[0];
			if(debug)
				System.out.println("Species: "+species+" conc: "+k1);
			if(context.symbols.contains(species)){
				if(debug)
					System.out.println("Done");
				init[context.symbols.indexOf(species)] = k1;
			}
			break;
		case KPLUS:
		case KMINUS:
		default:
			System.err.println("Unsupported parameter type: "+Constants.typeCurrentTarget1);
		}
		
		params = Constants.currentTarget2.split(" ");
		
		switch(Constants.typeCurrentTarget2){
		case KM:
			enz = params[0];
			if(params.length>=2){
				species = params[1];
			}
			trueEnz = context.enzymes.get(enz);
			if(debug)
				System.out.println("Enz: "+enz+" species: "+species+" trueEnz: "+trueEnz);
			if(trueEnz != null && (species.equals("") || trueEnz.Kms.get(species)!=null)){
				if(debug)
					System.out.println("Done");
				if(params.length>=2){
					trueEnz.Kms.put(species, k2);
				} else { //we actually update everyone, with a multiplying factor
					for(String s : trueEnz.Kms.keySet()){
						//find the original value we stored
						double defaultValue = context.Kmvalues.get(context.enzymenames.indexOf(enz)).get(context.symbolList.indexOf(s));
						//multiply by the rangeValue
						trueEnz.Kms.put(s, k2*defaultValue);
					}
				}
			}
			break;
		case KI:
			enz = params[0];
			if(params.length>=2){
				species = params[1];
			}
			trueEnz = context.enzymes.get(enz);
			if(debug)
				System.out.println("Enz: "+enz+" species: "+species+" trueEnz: "+trueEnz);
			if(trueEnz != null && (species.equals("") || trueEnz.kis.get(species)!=null)){
				if(debug)
					System.out.println("Done");
				if(params.length>=2){
					trueEnz.kis.put(species, k2);
				} else { //we actually update everyone, with a multiplying factor
					for(String s : trueEnz.Kms.keySet()){
						//find the original value we stored
						double defaultValue = context.Kmvalues.get(context.enzymenames.indexOf(enz)).get(context.symbolList.indexOf(s));
						//multiply by the rangeValue
						trueEnz.kis.put(s, k2*defaultValue);
					}
				}
			}
			break;
		case CONC:
			species = params[0];
			if(debug)
				System.out.println("Species: "+species+" conc: "+k2);
			if(context.symbols.contains(species)){
				if(debug)
					System.out.println("Done");
				init[context.symbols.indexOf(species)] = k2;
			}
			break;
		case KPLUS:
		case KMINUS:
		default:
			System.err.println("Unsupported parameter type: "+Constants.typeCurrentTarget1);
		}
		
	}

/*
	static public void displayBifurcationDiagram ( Individual_CR individual )	
	{
	
		Individual_CR indClone = (Individual_CR) individual.clone(); // build a clone that we can toy with.
		
		// [n]: IMPORTANT NOTE: I (nicolas) assume that only template initial concentration, K_m (and possibly k_i) can be modified. k+ and k- are *not* modified.
		
		//Lexing the equations
		double k1 = 0.1;
		double maxk1 = 4.0;
		int sizex = 5;
		int sizey = 5;
		double mink4 = 0.0;
		double maxk4 = 8.0;
		double stepx = (maxk4-mink4)/((double)sizex);
		double stepy = (maxk1-k1)/((double)sizey);
	
		String exp = "";
		JFrame frame = new JFrame();
		JPanel content = new JPanel(new GridLayout(0,sizex));
		double time = System.currentTimeMillis();
		
		double[][] matrixKm = indClone.getCopyOfMatrixKm(); // get Km values 
		double[][] matrixKi = indClone.getCopyOfMatrixKi(); // get Ki values
		
		while(k1 < maxk1){
			double k4 = mink4;
			while(k4 < maxk4){
				
				matrixKm[1][0] = 1.0/(k1*8.7e-2); // overwrite
				matrixKi[1][0] = 1.0/8.7e-2; // overwrite
				
				indClone.setMatrixKm(matrixKm);
				indClone.setMatrixKi(matrixKi);
				
				crnLexer myLex = new crnLexer(new ANTLRInputStream( indClone.getDescription() ) );
	
				//Parsing the equations
				crnParser myPars = new crnParser(new CommonTokenStream(myLex));
				//Parsing result
				//Contains "symbols": X,Y,Z... list of all substrats.
				//		   "enzymes": E1... list of all enzymes.
				//         "reactions": list of all reactions.
				ReactionsetContext context = myPars.reactionset();
	
				//initial concentrations for substrats: 10 nM all around!
				double[] init = new double[context.symbols.size()];
				Arrays.fill(init, 5);
	
				double[][] result = solveEqus(init,context.symbols,context.enzymes,context.reactions);
				//System.out.println(k1+" "+k4);
				fitness.OscillatorBasedResult fit = ODESineFitness.evaluateFitness(result, context.symbols);
				exp += (fit.finalResult()>1?" "+fit.average+"\t":" -\t");
				String[] seqNames = new String[context.symbols.size()];
				seqNames =context.symbols.toArray(seqNames);
	
				PlotExpData demo = new PlotExpData("",
						result, seqNames, "", "");
	
				content.add(demo);
	
				k4 += stepx;
			}
			//System.out.println(exp);
			exp += "\n";
			k1+= stepy;
		}
		double time2 = System.currentTimeMillis();
		System.out.println(exp);
		System.out.println("\n"+(time2-time));
		frame.add(content);
		frame.pack();
		frame.setVisible(true);
	}
	*/
	
	static public void displayBifurcationDiagram ( String s) 
	{
		boolean debug = true;
		boolean debug_verbose = false; //display the parameter settings
		
		int sizex = 5;
		int sizey = 5;
		
		
		JFrame frame = new JFrame();
		JPanel content = new JPanel(new GridLayout(0,sizex));
		double time = System.currentTimeMillis();
		
		//Lexing the equations
		crnLexer myLex = new crnLexer(new ANTLRInputStream(s));
		//Parsing the equations
		crnParser myPars = new crnParser(new CommonTokenStream(myLex));
		//Parsing result
		//Contains "symbols": X,Y,Z... list of all substrats.
		//		   "enzymes": E1... list of all enzymes.
		//         "reactions": list of all reactions.
		ReactionsetContext context = myPars.reactionset();
		
		//initial concentrations for substrats: 10 nM all around!
		double[] init = new double[context.symbols.size()];
		Arrays.fill(init, 5);
		
		if ( debug )
			System.out.println(s);
		
		// [n]: IMPORTANT NOTE: I (nicolas) assume that only template initial concentration, K_m (and possibly k_i) can be modified. k+ and k- are *not* modified.
		
		
		double k1 = Constants.range1min;
		double k2 = Constants.range2min;
		double stepx = (Constants.range2max-Constants.range2min)/((double) sizex);
		double stepy = (Constants.range1max-Constants.range1min)/((double) sizey);
		
		// compute phase diagram
		
		double[][] bifurcationDiagram = new double[ (int) (( Constants.range2max - Constants.range2min ) / stepx) ][ (int)(( Constants.range1max - Constants.range1min ) / stepy) ];
		int i = 0;
		int j = 0;
			
		while(k1 < Constants.range1max){
			
			k2 = Constants.range2min;
			i=0;			
			
			while(k2 < Constants.range2max){
				
				//overwrite params
				Arrays.fill(init, 5);
				setParams(context, init, k1, k2);
				
				if(debug_verbose){
					System.out.println("Ki :"+context.enzymes.get("E1").kis.get("X")+" Y: "+context.enzymes.get("E1").kis.get("Y")+" Z: "+context.enzymes.get("E1").kis.get("Z"));
					System.out.println("Km :"+context.enzymes.get("E1").Kms.get("X")+" Y: "+context.enzymes.get("E1").Kms.get("Y")+" Z: "+context.enzymes.get("E1").Kms.get("Z"));
					System.out.println("Conc :"+context.symbols+" "+init[0]+" "+init[1]+" "+init[2]);
				}
				double[][] result = solveEqus(init,context.symbols,context.enzymes,context.reactions);
				
				fitness.OscillatorBasedResult fit = ODESineFitness.evaluateFitness(result, context.symbols);
				
				bifurcationDiagram[i][j] = fit.finalResult();//>1?fit.average:0;
				
				String[] seqNames = new String[context.symbols.size()];
				seqNames =context.symbols.toArray(seqNames);
	
				PlotExpData demo = new PlotExpData("",
						result, seqNames, "", "");
	
				content.add(demo);
				k2 += stepx;
				i++;
			}
			
			k1 += stepy;
			j++;
		}

		// compute fitness (compare phase diagram with target phase diagram

		// TODO: toy fitness, target is flat everywhere
		double fitnessValue = 0; //Double.MAX_VALUE;
		for ( j = 0 ; j != bifurcationDiagram[0].length ; j++ )
		{
			for ( i = 0 ; i != bifurcationDiagram.length ; i++ ){
				fitnessValue+= bifurcationDiagram[i][j];
			}
		}
		
		// debug
		if ( debug )
		{
			System.out.println("Bifurcation diagram: ");
			for ( j = 0 ; j != bifurcationDiagram[0].length ; j++ )
			{
				System.out.print(" ");
				for ( i = 0 ; i != bifurcationDiagram.length ; i++ )
					System.out.print(bifurcationDiagram[i][j]+"\t");
				System.out.println();
			}
			System.out.println("Fitness value: " + fitnessValue +"\n");
		}
		
		double time2 = System.currentTimeMillis();
		
		System.out.println("\n"+(time2-time));
		frame.add(content);
		frame.pack();
		frame.setVisible(true);
		
	}
// =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
	
	
	public static void main(String[] args){
		double time = System.currentTimeMillis();
		
		/*getFitnessValue_2(
		//displayBifurcationDiagram(		
				"enzymes: E1 E2 E3"+
				"products: X Y Z"+
				"K_m[3][3] : { { 0.6743507837577559 0.6358068375835202 1.268509077624075 }{ 0.8149509950676753 1.9114570550215748 0.21570359817855497 }{ 4.51958163518514 0.9325297467992038 0.1070592639280451 } }"+
				"k_i[3][3] : { { 0.961387267149697 0.2454697693566046 0.3784668045256365 }{ 0.10157912376611032 5.2290560395547665 0.9875062488193197 }{ 0.13170194315139835 4.064847951618307 0.5054319675122235 } }"+
				"E1 + X ==> E1 + X + X"+
				"X + Y ==> Y + Y [ k+ = 0.139367940620003 ; k- = 0.0 ]"+
				"E2 + Y ==> E2"+
				"Z ==>  [ k+ = 0.5742314986545423 ; k- = 0.0 ]"
			);*/
		
		displayBifurcationDiagram(		
				"enzymes: E1 E2 E3"+
				"products: X Y Z"+
				"K_m[3][3] : { { 0.18765339747995985 0.16875402202331313 0.1220542785050846 }{ 0.31102316617374043 1.2231403074237865 2.183775345955116 }{ 0.16370457405019204 0.8140063918011246 0.3867744824566277 } }"+
				"k_i[3][3] : { { 5.251425290346123 0.13707797225592214 7.9447719273687465 }{ 4.7111404524971645 0.4543646305645125 0.24858524193583262 }{ 0.29351121487714205 1.8162883597662227 0.1320026714178671 } }"+
				"E2 + X ==> E2"+
				"Z + Z <==  [ k+ = 0.0 ; k- = 1.2230295280555032 ]"+
				"Y + X ==> Z + X [ k+ = 0.5251076010311138 ; k- = 0.0 ]"+
				"X + X <==> Z + X [ k+ = 0.2278018163007697 ; k- = 3.3514925643037796 ]");

		//getFitnessValue_1( Global_CR.returnDefaultIndividual(null) );		
		//displayBifurcationDiagram( Global_CR.returnDefaultIndividual(null) );
		//testDefault();
		
		double time2 = System.currentTimeMillis();
		System.out.println("\nCompleted in "+(time2-time)+"ms");
		
	}
	
}
